setwd('~/Dropbox/Abortion/Pre_Roe_Stuff/Paper Files/JLC submission/Replication Files')
#GETTING R SETUP.
#HERE IS VERSION OF R I USED, ALONG WITH VERSIONS OF PACKAGES USED
# > sessionInfo()
# R version 3.1.0 (2014-04-10)
# Platform: x86_64-apple-darwin10.8.0 (64-bit)
# 
# locale:
# [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# 
# attached base packages:
# [1] grid      stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#  [1] apsrtable_0.8-8 mrp_1.0-1       mrpdata_1.0     reshape2_1.4    plyr_1.8.1      lattice_0.20-29
#  [7] sp_1.0-15       blme_1.0-2      foreign_0.8-61  arm_1.6-10      lme4_1.1-6      Rcpp_0.11.3    
# [13] Matrix_1.1-3    MASS_7.3-31    
# 
# loaded via a namespace (and not attached):
# [1] abind_1.4-0         coda_0.16-1         minqa_1.2.3         nlme_3.1-117       
# [5] RcppEigen_0.3.2.2.0 splines_3.1.0       stringr_0.6.2      

#LOAD LIBRARIES
library(arm)
library(foreign)
library(blme)
library(mrp)
library(mrpdata)
library(apsrtable)
###################################
#TO INSTALL MRP PACKAGE FROM GITHUB, RUN THESE COMMANDS
#library(devtools) ; install_github("mrp", "malecki", sub="mrpdata")
#library(devtools) ; install_github("mrp", "malecki", sub="mrp")

rm(list=ls(all=TRUE))

#function for standard error of a mean
mean.se <- function (x, na.rm) {stderr <- sd(x, na.rm = TRUE) 
	return(stderr/sqrt(sum(!is.na(x))))}

#THESE SCRIPTS CREATE THE FUNCTIONS TO ESTIMATE UNCERTAINTY FROM THE MRP ESTIMATES
	
source("poststrat_w_uncertainty.R")
source("fitted.R")

options(digits = 3)
options(scipen = 100)

#number of sims for uncertainty
n.sims <- 1000
lower.ci.bound <- .025
upper.ci.bound <- .975

#SET UP STATELEVEL For MRP
statelevel <- read.dta("statelevel_to_use_in_mrp.dta", convert.underscore=T)
statelevel$state <- as.factor(statelevel$sstate )
statelevel <- statelevel[order(statelevel$state),]

#just keep variables I use
statelevel <- subset(statelevel, select=c(state, sstatename, percent.catholic.mean, mcgovern72, state.abortion.index, abortion.policy.change
		, heard.case, strike.overall, strike.any, strike.overall, urban.percent, fed.case))
		
#now read in "Synthetic" census using Lucas' method
census.1970.synthetic <- read.dta("1970_poststrat_race_wb_synthetic.dta", convert.underscore=T)
census.1970.synthetic <- census.1970.synthetic[ order(census.1970.synthetic$cstate),]
census.1970.synthetic$crace.wb <- factor(ifelse(census.1970.synthetic$crace.wb =="1","white","black"))

## To ensure matching of strata between poll and population,
## both should be factors with identical names and levels.
census.1970.synthetic <- within(census.1970.synthetic,{
    age.cat <- factor(cage.cat)
		edu.cat <- factor(cedu.cat)
		female <- factor(cfemale, labels=c("Male","Female"), exclude=NA)
		state <- factor(cstate)
		region <- factor(cregion)
		race.wb <- factor(crace.wb)
		cath <- factor(ccath, labels=c("Not Cath", "Cath"))
		f.race <- interaction(female,race.wb)
		f.race.cath <- interaction(female,race.wb, cath)
		cf.race <- interaction(female,race.wb)
		crace.female <- ( as.numeric(cfemale) *3) + as.numeric(crace.wb)
		cf.race.cath <- interaction(female,race.wb, cath)
		cage.edu.cat <- 4 * ( as.numeric(cage.cat) -1) + as.numeric(cedu.cat)
		})
	
########################################################################################################
#1st) DESCRIBE POLICY
########################################################################################################
#barplot of policy change
wrap.it <- function(x, len)
{ 
  sapply(x, function(y) paste(strwrap(y, len), collapse = "\n"),  USE.NAMES = FALSE)
}
 #Call this function with a list or vector
wrap.labels <- function(x, len)
{
  if (is.list(x))
  {
    lapply(x, wrap.it, len)
  } else {
    wrap.it(x, len)
  }
}

foo <- barplot(table(statelevel$state.abortion.index), axes=F, xlab="", horiz=F, density=0)
dev.off()
axis.text <- 1.2
label.text <- 1.5
text.size <- .7

pdf("Figure_1.pdf", height=5.5, width=8)

par(mar=c(3.2,4,3,.5), mfrow=c(1,1))

barplot(table(statelevel$state.abortion.index), axes=F, xlab="", horiz=F, density=0, names.arg=F)
polygon (x=c(2.5,2.5, par()$usr[2], par()$usr[2]),y=par()$usr[c(3,4,4,3)], border="NA", col="gray70")
barplot(table(statelevel$state.abortion.index),  density=0,add=T,axes=F, names.arg=F)

axis(1,labels=c(0:5), at =foo, tcl=F, mgp=c(2,.3,0), cex.axis=axis.text)
axis(2, las=1, cex.axis=axis.text)

#add state names inside bars
simpleCap <- function(x) {
  s <- strsplit(x, " ")[[1]]
  paste(toupper(substring(s, 1,1)), substring(s, 2),
      sep="", collapse="")#adjust collapse for spacing between names
}
temp.names <- sapply(statelevel$sstatename,simpleCap)
for (j in 0:max(statelevel$state.abortion.index) ){
		ok <- statelevel$state.abortion.index==j
		t.names <-  paste(temp.names[ok], " ", collapse="")
		wr.lap <- wrap.labels(t.names, length(t.names))
		#adjust two-word states
		wr.lap <- ifelse(	wr.lap=="Louisiana\nNewHampshire\nPennsylvania", "Louisiana\nNew Hampshire\nPennsylvania",wr.lap)
 		wr.lap <- ifelse(	wr.lap=="Alaska\nHawaii\nNewYork\nWashington", "Alaska\nHawaii\nNew York\nWashington",wr.lap)
 		wr.lap <- ifelse(	wr.lap=="Alabama\nArizona\nConnecticut\nD.c.\nIowa\nIdaho\nIllinois\nIndiana\nKentucky\nMassachusetts\nMaine\nMichigan\nMinnesota\nMissouri\nMontana\nNorthDakota\nNebraska\nNewJersey\nNevada\nOhio\nOklahoma\nRhodeIsland\nSouthDakota\nTennessee\nTexas\nUtah\nVermont\nWisconsin\nWestVirginia\nWyoming", "Alabama\nArizona\nConnecticut\nIowa\nIdaho\nIllinois\nIndiana\nKentucky\nMassachusetts\nMaine\nMichigan\nMinnesota\nMissouri\nMontana\nNorth Dakota\nNebraska\nNew Jersey\nNevada\nOhio\nOklahoma\nRhode Island\nSouth Dakota\nTennessee\nTexas\nUtah\nVermont\nWashington DC\nWisconsin\nWest Virginia\nWyoming",wr.lap)
		wr.lap <- ifelse(	wr.lap=="California\nColorado\nDelaware\nKansas\nMaryland\nNorthCarolina\nNewMexico\nOregon\nSouthCarolina\nVirginia", "California\nColorado\nDelaware\nKansas\nMaryland\nNorth Carolina\nNew Mexico\nOregon\nSouth Carolina\nVirginia",wr.lap)
 	 	text(foo[j+1], 0-.3, wr.lap ,pos=3, cex=text.size)
		}
#add labels
text(1.2, 31, "No Policy Change",xpd=NA,font=1,cex=1.5)
text(4.98,31, "Policy Change",xpd=NA, font=1,cex=1.5)
mtext("Index of liberalization of state policy", 1, line=1.9, cex=label.text, font=2)
mtext("Number of states", 2, line=2.5, cex=label.text, font=2)
box(bty="l")
box(which="outer")
mtext("Abortion policy in the states", 3, line=1.4, cex=1.8, font=2)

dev.off()	

####
#AGGREGATE OPINION
###
agg.results <- read.dta("aggregate_index_results_1962_1972.dta")
#plot
x <- unique(agg.results$year)
title.text <- 2
axis.text <-1.2
line.label.text <- 1.2

#note: pooling no and don't knows at the moment.
pdf("Figure_2.pdf", height = 5, width=7)

par(mfrow=c(1,1), mar=c(2,3,4,3))

plot(x, agg.results$mother_should_be_legal, ylim=c(0,100), type="b", lty=1, 
	pch=19, axes=F, xaxs="r", yaxs="i", col="red")
axis(1,at=x, mgp=c(2,.5,0), cex.axis = axis.text)
axis(2, at =c(0,25,50,75,100), las=1, mgp=c(2,.5,0), cex.axis = axis.text)
axis(4, at =c(0,25,50,75,100), las=1, mgp=c(2,.5,0), cex.axis = axis.text)

points(x,  agg.results$deformed_should_be_legal, ylim=c(0,100), type="b", lty=2, col="black", pch=19)
points(x,  agg.results$money_should_be_legal, ylim=c(0,100), type="b", lty=3, col="blue", pch=19)
points(x,  agg.results$enoughchildren_should_be_legal, ylim=c(0,100), type="b", lty=4, col="dark green", pch=19)
mtext("Abortion views, 1962-1972", 3, line=2, cex=title.text,font=2)
mtext("Percent saying abortion should be legal...", line=.5, cex=1.5, font=3)
text(1966, 82, "if the health of the mother is in danger",col="red",cex=line.label.text)
text(1966, 63, "where the child may be born deformed",  srt=5,cex=line.label.text)
text(1962, 27, "where the family does not have enough\nmoney to support another child", col="blue", adj=0,cex=line.label.text)
text(1972, 27, "where the parents\nsimply have all the\nchildren they want", col="dark green", adj=1, srt=30,cex=line.label.text)
box()

dev.off()

##########
#BEGIN  WORKING WITH INDIVIDUAL_LEVEL OPINION
########
#READ in MEGAPOLL
megapoll <- read.dta("pre_roe_megapoll.dta", convert.underscore=T, convert.factors = F)
dim(megapoll)
megapoll$state <- ifelse(megapoll$state=="",NA, megapoll$state)
megapoll <- megapoll[order(megapoll$state),]
table(megapoll$state) #MISSING AK, HI, ND, NM, and NV
#drop observations with missing observations
megapoll <- subset(megapoll, !is.na(state) & !is.na(age.cat) & !is.na(edu.cat) & !is.na(female) & !is.na(female)
	& !is.na(region) * !is.na(catholic))
megapoll$sstate <- NULL
megapoll$sstate.initnum <- NULL
#DROP 1962 poll
megapoll <- subset(megapoll, poll!="gallup_1962_08_23_index")

## To ensure matching of strata between poll and population,
## both should be factors with identical names and levels.
megapoll  <- within (megapoll , {
  age.cat <- factor(age.cat,exclude=NA)
  edu.cat <- factor(edu.cat,exclude=NA)
  female <- factor(female, labels=c("Male","Female"), exclude=NA)
  race.wb <- factor(race.wb,labels=c("white","black"),exclude=NA) 
  state <- factor(state,exclude=NA) 
  region <- factor(region,exclude=NA) 
  f.race <- interaction(female,race.wb)
  cath <- factor(catholic,  labels=c("Not Cath", "Cath"), exclude=NA)
  f.race.cath <- interaction(female,race.wb, cath)
})

#SEPARATE FOR TWO ANALYSES
#NB: I use "early" to denote the full repal measure, meaning you can get 
	#an abortion without restrictions in first 3 months (i.e. "early")
megapoll.for.SQ <- subset(megapoll, poll.type=="status quo" & !is.na(change.SQ))#remove all missing data from DV
megapoll.for.SQ$state <- droplevels(megapoll.for.SQ$state)#drop missing levels
megapoll.for.early <- subset(megapoll, poll.type=="3months" & !is.na(allow.early)) #remove all missing data from DV
mean(megapoll.for.early$allow.early,na.rm=T)
megapoll.for.early$state <- droplevels(megapoll.for.early$state)#drop missing levels
	#half the states have les than 40 respondens
#FOOTNOTE 15: check monotonicy of responses to status quo questions
responses <- cbind(megapoll.for.SQ$abort.mother,megapoll.for.SQ$abort.deformed,
									megapoll.for.SQ$abort.money,megapoll.for.SQ$abort.enough.children)
								
string <-  paste( responses[,1] , responses[,2], responses[3,], responses[,4])

table(responses[,1] - responses[,2] )#so all errors are people who do not support abortion to 
																		 #save mother but do where child may be deformed
table(responses[,1] - responses[,3] )
table(responses[,1] - responses[,4] )
table(responses[,2] - responses[,3] )
table(responses[,2] - responses[,4] )
table(responses[,3] - responses[,4] )
#More than 90% of responses exhibited transitive preferences.
non.mon <- ifelse(string=="0 0 1 0" | string =="0 1 0 0" | string =="0 1 0 NA" | string=="0 1 1 0" | string=="0 1 NA 0" | string=="0 NA 1 0" 
	|	string=="0 NA 1 NA" | string=="1 0 0 1" | string=="1 0 1 0" | string==" 1 0 1 1" | string=="1 0 1 NA" | string=="1 1 0 1" | string=="1 NA 0 1" 
	| string=="NA 0 1 0" | string=="NA 0 1 NA" | string==" NA 0 NA 1",1,0)
mean(non.mon, na.rm=T)	
megapoll.for.SQ.non.mom <- subset(megapoll.for.SQ ,non.mon!=1 )
#BELOW I check correlation when you drop non.monotonic responses 

#########################################################################
#individual-level models
#PRESENTED in Table A-1 of the appendix
logit.individual.sq <- glm(change.SQ ~
		as.factor(age.cat) + as.factor(edu.cat) +	female +	race.wb		+ catholic	+ dem + gop	+ as.factor(poll)
		,data = megapoll.for.SQ, family=binomial(link="logit"))
summary(logit.individual.sq)
logit.individual.early <- glm(allow.early ~
		as.factor(age.cat) + as.factor(edu.cat) +	female +	race.wb		+ catholic	+ dem + gop	+ as.factor(poll)
		,data = megapoll.for.early, family=binomial(link="logit"))
summary(logit.individual.early, detail=T,digits=2)

coefs <- c("Intercept", "Age 30-44", "Age 45-64", "Age 65+", 
		"Grade 9-11", "HS Grad", "Some college+", "Female", "Black", "Catholic",
		"Democrat", "Republican", "Gallup 1969", "NES", "Gallup 1972")
ind.reg.table <- apsrtable(logit.individual.sq,  logit.individual.early
	,	coef.names = coefs, model.names=c("Status quo", "Full repeal") , label="ind_reg_table"
	, caption.position="below"	, digits=2	, caption="Logit models of individual-level opinion. The first model looks at opinion on changing the statue quo, the second model looks at opinion on full repeal."
	)
ind.reg.table

library(DAMisc) #CALCUATE PRE
c(mean(logit.individual.sq$y), mean(logit.individual.early$y))
c(pre(logit.individual.sq)$pcp , pre(logit.individual.early)$pcp ) 
c(pre(logit.individual.sq)$pre, pre(logit.individual.early)$pre)
detach("package:DAMisc", unload=TRUE)                                                                
################
#BEGIN STATE-LEVEL ESTIMATES
##########################

#CALCUATE DISAGRREGATED PREDICTION FOR SQ MESAURE
sq.disagg.pred <- tapply(megapoll$change.SQ,megapoll$state, mean, na.rm=T)
sq.disagg.pred <- sq.disagg.pred[!is.nan(sq.disagg.pred)]
#UNCERTAINTY
sq.disagg.se <- tapply(megapoll$change.SQ,megapoll$state, mean.se, na.rm=T)[!is.nan(sq.disagg.pred)]
sq.disagg.se <- sq.disagg.se[!is.na(sq.disagg.se)]
#sim using rnorm()
sq.disagg.sim <- matrix(NA, nrow=length(sq.disagg.pred), ncol=n.sims)
rownames(sq.disagg.sim ) <- names(sq.disagg.pred)
for (i in 1:length(sq.disagg.pred)){
	sq.disagg.sim[i,] <- rnorm(n.sims,sq.disagg.pred[i], sq.disagg.se[i])
}
sq.disagg.sim <- ifelse(sq.disagg.sim<0,0,sq.disagg.sim)#truncate values over 0 and 1
sq.disagg.sim <- ifelse(sq.disagg.sim>1,1,sq.disagg.sim)

#SQ MRP MODEL
sq.mrp <- mrp(change.SQ ~ state+age.cat+edu.cat + f.race.cath
                  , data=megapoll.for.SQ,
                  population=census.1970.synthetic,
                  pop.weights="cpercent.state",
 				  				formula.pop.update= .~.-poll,
                  grouplevel.data.frames=list(statelevel),
    			  			formula.model.update=  .~.  + urban.percent
									)
											
sq.model <- getModel(sq.mrp)
summary(sq.model, detail=T, digits=4)
coef(sq.model)
se.coef(sq.model)
fixef    (sq.model)
se.fixef(sq.model)
sq.mrp.pred <- poststratify(sq.mrp , ~ state)


######################################################
#check correlation when you drop non.monotonic responses 
sq.mrp.non.mon <- mrp(change.SQ ~ state+age.cat+edu.cat + f.race.cath, data=megapoll.for.SQ.non.mom, population=census.1970.synthetic,
                  pop.weights="cpercent.state",	formula.pop.update= .~.-poll,  
									grouplevel.data.frames=list(statelevel),formula.model.update= .~.  + urban.percent)				
sq.mrp.pred.non.mon <- poststratify(sq.mrp.non.mon  , ~ state)
cor(sq.mrp.pred,sq.mrp.pred.non.mon )#.99
#check correlation when you drop NES
megapoll.for.SQ.drop.NES <- subset(megapoll, poll.type=="status quo" & !is.na(change.SQ) & poll!="NES")#remove all missing data from DV
sq.mrp.drop.NES <- mrp(change.SQ ~ state+age.cat+edu.cat + f.race.cath, data=megapoll.for.SQ.drop.NES, population=census.1970.synthetic,
                  pop.weights="cpercent.state",	formula.pop.update= .~.-poll,  
									grouplevel.data.frames=list(statelevel),formula.model.update= .~.  + urban.percent)				
sq.mrp.pred.drop.NES <- poststratify(sq.mrp.drop.NES  , ~ state)
cor(sq.mrp.pred,sq.mrp.pred.drop.NES ) #(.93)

#compare to disagg means
temp1 <- data.frame(state=names(sq.disagg.pred), sq.disagg.pred=as.numeric(sq.disagg.pred))
temp2 <- data.frame(state=names(sq.mrp.pred  ), sq.mrp.pred=as.numeric(sq.mrp.pred ))
sq.compare <- join(temp1, temp2)
sq.compare <- subset(sq.compare, !is.na(sq.disagg.pred))
cor(sq.compare$sq.disagg.pred, sq.compare$sq.mrp.pred) 
#calculate uncertainty
sq.mrp.sim <- poststratify.sim(sq.mrp  , ~  state,n.sims)  
sq.mrp.sim.quantiles <- apply(sq.mrp.sim ,1, quantile, probs=c(lower.ci.bound, .5,upper.ci.bound))#create matrix with medians and confidence intervals
#standard deviation
sq.mrp.sim.sd <- apply(sq.mrp.sim ,1, sd)

#plot of disagg vs. predicted among those states (Figure A-1 in the Appendix)
pdf("Figure_A1_A.pdf", height=5,width=6)

axis.text <- 1.2
title.text <-1.4
par(mar=c(3,3,1,1))
plot(sq.compare$sq.disagg.pred, sq.compare$sq.mrp.pred, type="n",  xlim=c(0,1), ylim=c(0,1), xaxs="i", yaxs="i", axes=F,xlab="",ylab="")
text(sq.compare$sq.disagg.pred, sq.compare$sq.mrp.pred, sq.compare$state)
lines(lowess(  sq.compare$sq.disagg.pred, sq.compare$sqmrp.pred))
axis(1, mgp=c(2,.5,0), at=seq(0,1,.10), labels=seq(0,100,10))
axis(2, mgp=c(2,.5,0), at=seq(0,1,.10), labels=seq(0,100,10), las=1)
mtext("Disaggregation estimates of support for changing SQ",1, line=2, cex=axis.text)
mtext("MRP estimates of support for changing SQ",2, line=2, cex=axis.text)
mtext("A) Support for changing the status quo ",3, line=-1, cex=title.text)
abline(0,1, lty=2)
corr <- round(cor(sq.compare$sq.disagg.pred, sq.compare$sq.mrp.pred),2)
text(.8,.2, expression(paste(rho," = ")), adj=0, cex=1.4)
text(.87,.21, corr, adj=0,, cex=1.4)

dev.off()


####################################################
###NOW USE GALLUP QUESTIONS ON LEGALIZATION in 1st 3 months
####################################################
# TURNING TO A COMPLETELY DIFFERENT TOPIC...WOULD YOU FAVOR
# OR OPPOSE A LAW WHICH WOULD PERMIT A WOMAN TO GO TO A
# DOCTOR TO END PREGNANCY AT ANY TIME DURING THE FIRST THREE
# MONTHS*
early.disagg.pred <- tapply(megapoll.for.early$allow.early,megapoll.for.early$state, mean, na.rm=T)
early.disagg.pred <- early.disagg.pred[!is.nan(early.disagg.pred)]
early.disagg.se <- tapply(megapoll.for.early$allow.early,megapoll.for.early$state, mean.se, na.rm=T)[!is.nan(early.disagg.pred)]
early.disagg.se <- early.disagg.se[!is.na(early.disagg.se)]

#what is correlation btw sq and early disagg
disagg.data <- join( data.frame(state=names(early.disagg.pred), early.disagg.pred=as.numeric(early.disagg.pred)),
							data.frame(state=names(sq.disagg.pred), sq.disagg.pred=as.numeric(sq.disagg.pred)))
cor(disagg.data$early.disagg.pred, disagg.data$sq.disagg.pred,use="complete.obs")#only .35

early.disagg.pred <- tapply(megapoll.for.early$allow.early,megapoll.for.early$state, mean, na.rm=T)
early.disagg.pred <- early.disagg.pred[!is.nan(early.disagg.pred)]
early.disagg.se <- tapply(megapoll.for.early$allow.early,megapoll.for.early$state, mean.se, na.rm=T)[!is.nan(early.disagg.pred)]
early.disagg.se <- early.disagg.se[!is.na(early.disagg.se)]
#sim using rnorm()
early.disagg.sim <- matrix(NA, nrow=length(early.disagg.pred), ncol=n.sims)
rownames(early.disagg.sim ) <- names(early.disagg.pred)
for (i in 1:length(early.disagg.pred)){
	early.disagg.sim[i,] <- rnorm(n.sims,early.disagg.pred[i], early.disagg.se[i])
}
early.disagg.sim <- ifelse(early.disagg.sim<0,0,early.disagg.sim)#truncate values over 0 and 1
early.disagg.sim <- ifelse(early.disagg.sim>1,1,early.disagg.sim)


early.mrp <- mrp(allow.early ~ state+age.cat+edu.cat + f.race.cath
                  , data=megapoll.for.early,
                  population=census.1970.synthetic,
                  pop.weights="cpercent.state",
 				  				formula.pop.update= .~.-poll,
                  grouplevel.data.frames=list(statelevel),
    			  			formula.model.update= .~.  + urban.percent
										)

early.model <- getModel(early.mrp)
summary(early.model, detail=T, digits=4)
coef(early.model)
se.coef(early.model)
fixef    (early.model)
se.fixef(early.model)
early.mrp.pred <- poststratify(early.mrp , ~ state)
early.mrp.pred
cor(early.mrp.pred, sq.mrp.pred)

#compare to disagg means
temp1 <- data.frame(state=names(early.disagg.pred), early.disagg.pred=as.numeric(early.disagg.pred))
temp2 <- data.frame(state=names(early.mrp.pred  ), early.mrp.pred=as.numeric(early.mrp.pred ))
early.compare <- join(temp1, temp2)
early.compare <- subset(early.compare, !is.na(early.disagg.pred))
cor(early.compare$early.disagg.pred, early.compare$early.mrp.pred) 

#calculate uncertainty
early.mrp.sim <- poststratify.sim(early.mrp  , ~  state,n.sims) 
early.mrp.sim.quantiles <- apply(early.mrp.sim ,1, quantile, probs=c(lower.ci.bound, .5,upper.ci.bound))#create matrix with medians and confidence intervals
#standard deviation
early.mrp.sim.sd <- apply(early.mrp.sim ,1, sd)

#plot of disagg vs. predicted among those states 
pdf("Figure_A1_B.pdf", height=5,width=6)

axis.text <- 1.2
title.text <-1.4
par(mar=c(3,3,1,1), mfrow=c(1,1))
plot(early.compare$early.disagg.pred, early.compare$early.mrp.pred, type="n",  xlim=c(0,1), ylim=c(0,1), xaxs="i", yaxs="i", axes=F,xlab="",ylab="")
text(early.compare$early.disagg.pred, early.compare$early.mrp.pred, early.compare$state)
axis(1, mgp=c(2,.5,0), at=seq(0,1,.10), labels=seq(0,100,10))
axis(2, mgp=c(2,.5,0), at=seq(0,1,.10), labels=seq(0,100,10), las=1)
mtext("B) Support for full repeal",3, line=-1, cex=title.text)
mtext("Disaggregation estimates of support for full repeal",1, line=2, cex=axis.text)
mtext("MRP estimates of support for full repeal",2, line=2, cex=axis.text)
abline(lm( early.compare$early.mrp.pred~ early.compare$early.disagg.pred))
corr <- round(cor(early.compare$early.disagg.pred, early.compare$early.mrp.pred),2)
text(.8,.2, expression(paste(rho," = ")), adj=0, cex=1.4)
text(.87,.21, corr, adj=0,, cex=1.4)

dev.off()

#plot early MRP vs SQ disagg and mrp

pdf("Figure_A2_B.pdf", height=5,width=6)

axis.text <- 1.2
title.text <-1.4
par(mar=c(3,3,1,1), mfrow=c(1,1))
plot(disagg.data$sq.disagg.pred, disagg.data$early.disagg.pred, type="n",  xlim=c(0,.80), ylim=c(0,.80), xaxs="i", yaxs="i", axes=F,xlab="",ylab="")
text(disagg.data$sq.disagg.pred, disagg.data$early.disagg.pred, disagg.data$state)
axis(1, mgp=c(2,.5,0), at=seq(0,1,.10), labels=seq(0,100,10))
axis(2, mgp=c(2,.5,0), at=seq(0,1,.10), labels=seq(0,100,10), las=1)
mtext("Support for changing SQ",1, line=2, cex=axis.text)
mtext("Support for full repeal",2, line=2, cex=axis.text)
mtext("B) Disaggregation estimates ",3, line=-1, cex=title.text)
abline(lm( disagg.data$early.disagg.pred~ disagg.data$sq.disagg.pred))
corr <- round(cor(disagg.data$sq.disagg.pred, disagg.data$early.disagg.pred, use="complete.obs"),2)
text(.6,.05, expression(paste(rho," = ")), adj=0, cex=1.4)
text(.67,.05, corr, adj=0,, cex=1.4)

dev.off()

pdf("Figure_A2_A.pdf", height=5,width=6)

temp <- join( data.frame(sq.mrp.pred, state=names(sq.mrp.pred) ) ,  data.frame( early.mrp.pred, state=names( early.mrp.pred)) )
axis.text <- 1.2
title.text <-1.4
par(mar=c(3,3,1,1), mfrow=c(1,1))
plot(temp$sq.mrp.pred, temp$early.mrp.pred, type="n",  xlim=c(0,.80), ylim=c(0,.80), xaxs="i", yaxs="i", axes=F,xlab="",ylab="")
text(temp$sq.mrp.pred, temp$early.mrp.pred, temp$state)
axis(1, mgp=c(2,.5,0), at=seq(0,1,.10), labels=seq(0,100,10))
axis(2, mgp=c(2,.5,0), at=seq(0,1,.10), labels=seq(0,100,10), las=1)
mtext("Support for changing SQ",1, line=2, cex=axis.text)
mtext("Support for full repeal",2, line=2, cex=axis.text)
mtext("A) MRP estimates ",3, line=-1, cex=title.text)
abline(lm( temp$early.mrp.pred~ temp$sq.mrp.pred))
corr <- round(cor(temp$sq.mrp.pred, temp$early.mrp.pred),2)
text(.6,.05, expression(paste(rho," = ")), adj=0, cex=1.4)
text(.67,.05, corr, adj=0,, cex=1.4)

dev.off()

#MAKE DOTPLOT OF SUPPORT, using MRP AND DISAGG
#SQ
sq.mrp.pred.sorted <- sort(sq.mrp.pred)
sq.mrp.pred.lower.ci <- sq.mrp.sim.quantiles[1,][order(sq.mrp.pred)]
sq.mrp.pred.upper.ci <- sq.mrp.sim.quantiles[3,][order(sq.mrp.pred)]
#EARLY
early.mrp.pred.sorted <- sort(early.mrp.pred)
early.mrp.pred.lower.ci <- early.mrp.sim.quantiles[1,][order(early.mrp.pred)]
early.mrp.pred.upper.ci <- early.mrp.sim.quantiles[3,][order(early.mrp.pred)]


###COMBINED DOTPLOTS
#MRP
pdf("Figure_3.pdf", height=7,width=10)

layout( rbind(c(1,1),c(2,3)), heights=c(.03,1))

par(mar=c(0,0,0,0))
plot(0,0, axes=F, type="n", main="", xlab="", ylab="")
text(0,-.6,"State-level public opinion on abortion policy", font=2, cex=1.8, xpd=NA)

y <- c(1:length(sq.mrp.pred ))
par(mar=c(3.8,2,.1,1))
plot(sq.mrp.pred.sorted, y, type="n", axes=F, xlim=c(.1,.9),pch=19, xlab="", ylab="", xaxs="i", yaxs="r")
axis(1, mgp=c(2,.5,0), at=seq(0,1,.10), labels=seq(0,100,10))
axis(2,at=y, label=names(sq.mrp.pred.sorted), las=1, cex.axis=.7, mgp=c(2,.5,0))
abline(v=50)
abline(h=y, lty=2, col="gray70")
segments(sq.mrp.pred.lower.ci,y, sq.mrp.pred.upper.ci,y)
points(sq.mrp.pred.sorted  , y, pch=19, bg="white")
mtext("Estimated level of support for changing status quo",1, line=2, cex=axis.text)
abline(v=.5, lty=2)

y <- c(1:length(early.mrp.pred ))
plot(early.mrp.pred.sorted, y, type="n", axes=F, xlim=c(.1,.9),pch=19, xlab="", ylab="", xaxs="i", yaxs="r")
axis(1, mgp=c(2,.5,0), at=seq(0,1,.10), labels=seq(0,100,10))
axis(2,at=y, label=names(early.mrp.pred.sorted), las=1, cex.axis=.7, mgp=c(2,.5,0))
abline(v=50)
abline(h=y, lty=2, col="gray70")
segments(early.mrp.pred.lower.ci,y, early.mrp.pred.upper.ci,y)
points(early.mrp.pred.sorted  , y, pch=19, bg="white")
mtext("Estimated level of support for full repeal",1, line=2, cex=axis.text)
abline(v=.5, lty=2)
		
dev.off()

##COMBINED DOTPLOTS:DISAGGREGRATION
#SQ
sq.disagg.pred.sorted <- sort(sq.disagg.pred)
sq.disagg.pred.lower.ci <- sq.disagg.pred-1.96*sq.disagg.se 
sq.disagg.pred.lower.ci <- sq.disagg.pred.lower.ci[order(sq.disagg.pred)]
sq.disagg.pred.upper.ci <- sq.disagg.pred+1.96*sq.disagg.se
sq.disagg.pred.upper.ci <- sq.disagg.pred.upper.ci[order(sq.disagg.pred)]
#EARLY
early.disagg.pred.sorted <- sort(early.disagg.pred)
early.disagg.pred.lower.ci <- early.disagg.pred-1.96*early.disagg.se 
early.disagg.pred.lower.ci <- early.disagg.pred.lower.ci[order(early.disagg.pred)]
early.disagg.pred.upper.ci <- early.disagg.pred+1.96*early.disagg.se
early.disagg.pred.upper.ci <- early.disagg.pred.upper.ci[order(early.disagg.pred)]

pdf("Figure_B1.pdf", height=7,width=10)

layout( rbind(c(1,1),c(2,3)), heights=c(.03,1))

par(mar=c(0,0,0,0))
plot(0,0, axes=F, type="n", main="", xlab="", ylab="")
text(0,-.6,"State-level public opinion on abortion policy", font=2, cex=1.8, xpd=NA)

y <- c(1:length(sq.disagg.pred ))
par(mar=c(3.8,2,.1,1))
plot(sq.disagg.pred.sorted, y, type="n", axes=F, xlim=c(0,1),pch=19, xlab="", ylab="", xaxs="i", yaxs="r")
axis(1, mgp=c(2,.5,0), at=seq(0,1,.10), labels=seq(0,100,10))
axis(2,at=y, label=names(sq.disagg.pred.sorted), las=1, cex.axis=.7, mgp=c(2,.5,0))
abline(v=50)
abline(h=y, lty=2, col="gray70")
segments(sq.disagg.pred.lower.ci,y, sq.disagg.pred.upper.ci,y)
points(sq.disagg.pred.sorted  , y, pch=19, bg="white")
mtext("Estimated level of support for changing status quo",1, line=2, cex=axis.text)
abline(v=.5, lty=2)

y <- c(1:length(early.disagg.pred ))
plot(early.disagg.pred.sorted, y, type="n", axes=F, xlim=c(0,1),pch=19, xlab="", ylab="", xaxs="i", yaxs="r")
axis(1, mgp=c(2,.5,0), at=seq(0,1,.10), labels=seq(0,100,10))
axis(2,at=y, label=names(early.disagg.pred.sorted), las=1, cex.axis=.7, mgp=c(2,.5,0))
abline(v=50)
abline(h=y, lty=2, col="gray70")
segments(early.disagg.pred.lower.ci,y, early.disagg.pred.upper.ci,y)
points(early.disagg.pred.sorted  , y, pch=19, bg="white")
mtext("Estimated level of support for full repeal",1, line=2, cex=axis.text)
abline(v=.5, lty=2)
		
dev.off()


######################################################################
#CONNNECT OPINION TO POLICY
######################################################################
#MERGE PUBLIC OPINION AND POLICY ONTO TO STATELEVEL DATA
#NEED TO MAKE SURE STATES ARE IN RIGHT ORDER
#give them different names 

#SO: MEASURES TO USE FOR LOOKING AT POLICY AND COURT DECISIONS
#any change: 
#favor legalization in 3 months:gallup.combined.3months.pred.lib 

#FOR ALL REGRESSIONS, USE UNCERTAINTY. MERGE SIMS UNTO STATELEVEL
#run regression with simulated errors
#you have n.sim ddisaggs of X
#for each sim, run regression.
#then sim those results and take one sim from the simulation of the regressions,
#Do n.sims time
temp.df <- data.frame(sq.mrp.sim)
temp.df$state <- rownames(temp.df)
statelevel.mrp <- join(statelevel,temp.df)
temp.df <- data.frame(sq.disagg.sim)
temp.df$state <- rownames(temp.df)
statelevel.disagg <- join(statelevel,temp.df)
#MRP
temp.df <- data.frame(sq.mrp.sim)
temp.df$state <- rownames(temp.df)
mrp.preds <- data.frame(state=names(sq.mrp.pred), early.mrp.pred=early.mrp.pred, sq.mrp.pred=sq.mrp.pred)
statelevel.mrp <- join(statelevel, mrp.preds)
statelevel.mrp <- join(statelevel.mrp, temp.df)
#DISAGG
temp.df <- data.frame(sq.disagg.sim)
temp.df$state <- rownames(temp.df)
statelevel.disagg <- join(statelevel, disagg.data)
statelevel.disagg <- join(statelevel.disagg, temp.df)

n.uncertain.lines <- 100
#create indicator for any policy change
statelevel.mrp$any.policy.change <- ifelse(statelevel.mrp$state.abortion.index >=2,1,0)
statelevel.disagg$any.policy.change <- ifelse(statelevel.disagg$state.abortion.index >=2,1,0)

############################
#Figure_4.pdf-- MRP: Policy vs. Opinion on Status Quo
pdf("Figure_4.pdf", height=4, width=8)

axis.text <- 1.2
title.text <- 1.4
y <- statelevel.mrp$state.abortion.index
x <- statelevel.mrp$sq.mrp.pred
summary(lm(y~x))

state.label.size <- .7
set.seed(3)
par(mfrow=c(1,1), mar=c(4,3,2,1))
#all states
ok <-  T
plot(x[ok],y[ok], type="n", ylim=c(0,5), xlim=c(.2,.8), axes=F, xlab="", ylab="")
axis(2, at = c(0,1,2,3,4,5), mgp=c(2,.5,0), cex.axis=axis.text)
axis(1,  at=seq(0,1,.10), labels=seq(0,100,10), mgp=c(2,.5,0), cex.axis=axis.text, las=1)
mtext("Estimated level of support for changing status quo",1, srt=90, line=2, cex=axis.text)
mtext("Liberalism of abortion policy",2, line=2, cex=axis.text)

#REGRESSION LINES WITH UNCERTAINTY
index.coefs <- matrix(NA, nrow=n.sims, ncol=2)#do intercept, then slope
for (i in 1:n.sims){
 temp.x <- which(names(statelevel.mrp)=="X1")
 foo <-	lm( statelevel.mrp$state.abortion.index ~ statelevel.mrp[,temp.x+i-1])#LOOP OVER COLUMNNS WITH SIMMED OPINION
 foo2 <- sim(foo, n.sims=1)
 index.coefs[i,1] <- coef(foo2)[1]#intercept
 index.coefs[i,2] <- coef(foo2)[2]#slope
}
quantile(index.coefs[,2],c(.025,.5,.975))
mean(index.coefs[,2]>0)

#get 95% confidence intervals at every level of opinion
opinion.sequence <- seq(.2,.8, .05)
index.preds <- matrix(NA, nrow=length(opinion.sequence), ncol=3)
for (j in 1:length(opinion.sequence)){
	index.preds[j,1] <- quantile(index.coefs[,1] + index.coefs[,2]*opinion.sequence[j], .025)
	index.preds[j,2]  <- quantile(index.coefs[,1] + index.coefs[,2]*opinion.sequence[j], .5)
	index.preds[j,3]  <- quantile(index.coefs[,1] + index.coefs[,2]*opinion.sequence[j], .975)
}

polygon(c(rev(opinion.sequence), opinion.sequence), c(rev(index.preds[ ,3]), index.preds[ ,1]), col = "gray70", border = NA)
points(opinion.sequence, index.preds[ ,2], col="black", type="l")
#add shading for change and no change
#polygon (x= c(par()$usr[2], par()$usr[2], par()$usr[1], par()$usr[1]) , y=c(0,1.5,1.5,0), border="NA", col="gray70")
mtext("Public opinion on changing status quo and abortion policy",3, line=.5, cex=title.text, font=2)
#first, states with change
ok <- statelevel.mrp$abortion.policy.change==1
#adjust jitter depending on policy
ok2 <- y==2 | y==4 
jitter.x <-.02
jitter.y <- .1
text(		x[ok & ok2] +rnorm(length(x[ok & ok2]),0,1.8*jitter.x),
	y[ok & ok2] +rnorm(length(x[ok & ok2]),0,jitter.y),
	statelevel.mrp$state[ok & ok2], 	srt=0, cex=state.label.size, col="dark green")
ok2 <- y!=2 & y!=4 
jitter.x <-.02
jitter.y <- .04
text(		x[ok & ok2] +rnorm(length(x[ok & ok2]),0,1.8*jitter.x),
	y[ok & ok2] +rnorm(length(x[ok & ok2]),0,jitter.y),
	statelevel.mrp$state[ok & ok2], 	srt=0, cex=state.label.size, col="dark green")
#then, states without change (do 0 and 1's separately)
ok <- statelevel.mrp$abortion.policy.change==0
ok2 <- y==1
jitter.x <-.03
jitter.y <- .15
text(		x[ok & ok2] +rnorm(length(x[ok & ok2]),0,1.8*jitter.x),
	y[ok & ok2] +rnorm(length(x[ok & ok2]),0,jitter.y),
	statelevel.mrp$state[ok & ok2], 	srt=0, cex=state.label.size, col="red")
ok <-  T
ok2 <- y==0
jitter.x <-0
jitter.y <- 0
text(		x[ok & ok2] +rnorm(length(x[ok & ok2]),0,1.8*jitter.x),
	y[ok & ok2] +rnorm(length(x[ok & ok2]),0,jitter.y),
	statelevel.mrp$state[ok & ok2], 	srt=0, cex=state.label.size, col="red")
ok <-  T

box()

dev.off()

#For text, what does a shift in policy predict 
#moving from 30 to 70
index.preds[3,2]
index.preds[11,2]

##Figure_B2.pdf-- DISAGG: Policy vs. Opinion on Status Quo
pdf("Figure_B2.pdf", height=4, width=8)

axis.text <- 1.2
title.text <- 1.4
y <- statelevel.disagg$state.abortion.index
x <- statelevel.disagg$sq.disagg.pred
summary(lm(y~x))

state.label.size <- .7
set.seed(3)
par(mfrow=c(1,1), mar=c(4,3,2,1))
#all states
ok <-  T
plot(x[ok],y[ok], type="n", ylim=c(0,5), xlim=c(0,1), axes=F, xlab="", ylab="")
axis(2, at = c(0,1,2,3,4,5), mgp=c(2,.5,0), cex.axis=axis.text)
axis(1,  at=seq(0,1,.10), labels=seq(0,100,10), mgp=c(2,.5,0), cex.axis=axis.text, las=1)
mtext("Estimated level of support for changing status quo (using disaggregation)",1, srt=90, line=2, cex=axis.text)
mtext("Liberalism of abortion policy",2, line=2, cex=axis.text)

#REGRESSION LINES WITH UNCERTAINTY
index.coefs <- matrix(NA, nrow=n.sims, ncol=2)#do intercept, then slope
for (i in 1:n.sims){
 temp.x <- which(names(statelevel.disagg)=="X1")
 foo <-	lm(statelevel.disagg$state.abortion.index~statelevel.disagg[,temp.x+i-1])#LOOP OVER COLUMNNS WITH SIMMED OPINION
 foo2 <- sim(foo, n.sims=1)
 index.coefs[i,1] <- coef(foo2)[1]#intercept
 index.coefs[i,2] <- coef(foo2)[2]#slope
}
quantile(index.coefs[,2],c(.025,.5,.975))
mean(index.coefs[,2]>0)

#get 95% confidence intervals at every level of opinion
opinion.sequence <- seq(0,1, .05)
index.preds <- matrix(NA, nrow=length(opinion.sequence), ncol=3)
for (j in 1:length(opinion.sequence)){
	index.preds[j,1] <- quantile(index.coefs[,1] + index.coefs[,2]*opinion.sequence[j], .025)
	index.preds[j,2]  <- quantile(index.coefs[,1] + index.coefs[,2]*opinion.sequence[j], .5)
	index.preds[j,3]  <- quantile(index.coefs[,1] + index.coefs[,2]*opinion.sequence[j], .975)
}

polygon(c(rev(opinion.sequence), opinion.sequence), c(rev(index.preds[ ,3]), index.preds[ ,1]), col = "gray70", border = NA)
points(opinion.sequence, index.preds[ ,2], col="black", type="l")
#add shading for change and no change
#polygon (x= c(par()$usr[2], par()$usr[2], par()$usr[1], par()$usr[1]) , y=c(0,1.5,1.5,0), border="NA", col="gray70")
mtext("Public opinion on changing status quo and abortion policy",3, line=.5, cex=title.text, font=2)
#first, states with change
ok <- statelevel.disagg$abortion.policy.change==1
#adjust jitter depending on policy
ok2 <- y==2 | y==4 
jitter.x <-.02
jitter.y <- .1
text(		x[ok & ok2] +rnorm(length(x[ok & ok2]),0,1.8*jitter.x),
	y[ok & ok2] +rnorm(length(x[ok & ok2]),0,jitter.y),
	statelevel.disagg$state[ok & ok2], 	srt=0, cex=state.label.size, col="dark green")
ok2 <- y!=2 & y!=4 
jitter.x <-.02
jitter.y <- .04
text(		x[ok & ok2] +rnorm(length(x[ok & ok2]),0,1.8*jitter.x),
	y[ok & ok2] +rnorm(length(x[ok & ok2]),0,jitter.y),
	statelevel.disagg$state[ok & ok2], 	srt=0, cex=state.label.size, col="dark green")
#then, states without change (do 0 and 1's separately)
ok <- statelevel.disagg$abortion.policy.change==0
ok2 <- y==1
jitter.x <-.03
jitter.y <- .15
text(		x[ok & ok2] +rnorm(length(x[ok & ok2]),0,1.8*jitter.x),
	y[ok & ok2] +rnorm(length(x[ok & ok2]),0,jitter.y),
	statelevel.disagg$state[ok & ok2], 	srt=0, cex=state.label.size, col="red")
ok <-  T
ok2 <- y==0
jitter.x <-0
jitter.y <- 0
text(		x[ok & ok2] +rnorm(length(x[ok & ok2]),0,1.8*jitter.x),
	y[ok & ok2] +rnorm(length(x[ok & ok2]),0,jitter.y),
	statelevel.disagg$state[ok & ok2], 	srt=0, cex=state.label.size, col="red")
ok <-  T

box()

dev.off()

#DO CONGRUENCE ANALYSIS WITH UNCERTAINTY
#just keep state, any policy change and sims
#MRP
temp.end <- paste0("X",n.sims)
temp1 <- statelevel.mrp[,which(names(statelevel.mrp)=="X1"): which(names(statelevel.mrp)==paste0("X",n.sims)) ]
any.policy.change <- statelevel.mrp[,"any.policy.change"]
state <- statelevel.mrp[,"state"]
temp.df <- data.frame(cbind(state, any.policy.change ,temp1))

congruence.data.policy <- melt(temp.df, id.vars = c("state", "any.policy.change"), 
	value.name="sq.sims", variable.name="sim")
congruence.data.policy <- congruence.data.policy[order(congruence.data.policy$state),]

#now create congruence variable
congruence.data.policy$congruent <- ifelse( (congruence.data.policy$sq.sims >=.5 & congruence.data.policy$any.policy.change == 1 ) |
	 (congruence.data.policy$sq.sims <.5 & congruence.data.policy$any.policy.change == 0 ),1,0 )
statelevel.mrp$state.cong.prop.sq <- tapply(congruence.data.policy$congruent, congruence.data.policy$state, mean)

#FIGURE 5
pdf("Figure_5.pdf", height=8, width=8)

x.axis.text <- 1.1
y.axis.text <- .7
label.text <- 1.4
title.text <- 1.6
state.text <- .8

prop.cong.no.change.sorted <- sort(statelevel.mrp$state.cong.prop.sq[statelevel.mrp$any.policy.change==0 ])

layout(c(1,2), heights=c(.8,1))
par(mar=c(4,1,1,1))

prop.cong.yes.change.sorted <- sort(statelevel.mrp$state.cong.prop.sq[statelevel.mrp$any.policy.change==1 ])
y <- c(length(prop.cong.yes.change.sorted):1)
plot(prop.cong.yes.change.sorted, y, axes=F, xlab="",ylab="", main="", pch=19, type="n")
text(prop.cong.yes.change.sorted, y, names(prop.cong.yes.change.sorted), cex=state.text)
axis(1, at=seq(0,1,.25), mgp=c(2,.5,0), cex.axis=x.axis.text)
#axis(2, at=y, labels=names(prop.cong.yes.change.sorted), las=1,mgp=c(2,.5,0), cex.axis=y.axis.text)
mtext("Probability that policy was congruent with opinion",1, line=2, cex=label.text)
mtext("States with policy change", 3, line=-2,adj=.9, cex=title.text,font=2)
box()

y <- c(length(prop.cong.no.change.sorted):1)
plot(prop.cong.no.change.sorted, y, axes=F, xlab="",ylab="", main="", pch=19, type="n")
text(prop.cong.no.change.sorted, y, names(prop.cong.no.change.sorted), cex=state.text)
axis(1, at=seq(0,1,.25), mgp=c(2,.5,0), cex.axis=x.axis.text)
#axis(2, at=y, labels=names(prop.cong.no.change.sorted), las=1,mgp=c(2,.5,0), cex.axis=y.axis.text)
mtext("Probability that policy was congruent with opinion",1, line=2, cex=label.text)
mtext("States with no policy change", 3, line=-3,cex=title.text,font=2)
box()

dev.off()

#CREATE 2x2 Table at Bottom of Figure 5
#another way -- do 2*2 table with uncertainty (as in KLMP). In each sim, get each cell. then build vec
no.policy.change.mrp.opinion.no.policy.change.mrp <- rep(NA, n.sims)
yes.policy.change.mrp.opinion.no.policy.change.mrp <- rep(NA, n.sims)
no.policy.change.mrp.opinion.yes.policy.change.mrp <- rep(NA, n.sims)
yes.policy.change.mrp.opinion.yes.policy.change.mrp <- rep(NA, n.sims)

for (i in 1:n.sims){
	ok <- congruence.data.policy$sim ==paste0("X", i)
 	foo <-	table( congruence.data.policy$any.policy.change[ok] , congruence.data.policy$sq.sims[ok] >= .5)
	no.policy.change.mrp.opinion.no.policy.change.mrp[i] <- foo[1] / length(congruence.data.policy$any.policy.change[ok])
	yes.policy.change.mrp.opinion.no.policy.change.mrp[i] <- foo[2] / length(congruence.data.policy$any.policy.change[ok])
	no.policy.change.mrp.opinion.yes.policy.change.mrp[i] <- foo[3] / length(congruence.data.policy$any.policy.change[ok])
	yes.policy.change.mrp.opinion.yes.policy.change.mrp[i] <- foo[4] / length(congruence.data.policy$any.policy.change[ok])
}
quantile(	no.policy.change.mrp.opinion.no.policy.change.mrp, c(.025,.5,.975))
quantile(	yes.policy.change.mrp.opinion.no.policy.change.mrp, c(.025,.5,.975))
quantile(	no.policy.change.mrp.opinion.yes.policy.change.mrp, c(.025,.5,.975))
quantile(	yes.policy.change.mrp.opinion.yes.policy.change.mrp, c(.025,.5,.975))


########
#SAME, BUT WITH DISAGG
temp.end <- paste0("X",n.sims)
#drop states without estimates
temp.statelevel.disagg <- subset(statelevel.disagg, !is.na(sq.disagg.pred))
temp1 <-temp.statelevel.disagg[,which(names(statelevel.disagg)=="X1"): which(names(statelevel.disagg)==paste0("X",n.sims)) ]
any.policy.change <- temp.statelevel.disagg[,"any.policy.change"]
state <- temp.statelevel.disagg[,"state"]
temp.df <- data.frame(cbind(state, any.policy.change ,temp1))

congruence.data.policy <- melt(temp.df, id.vars = c("state", "any.policy.change"), 
	value.name="sq.sims", variable.name="sim")
congruence.data.policy <- congruence.data.policy[order(congruence.data.policy$state),]

#now create congruence variable
congruence.data.policy$congruent <- ifelse( (congruence.data.policy$sq.sims >=.5 & congruence.data.policy$any.policy.change == 1 ) |
	 (congruence.data.policy$sq.sims <.5 & congruence.data.policy$any.policy.change == 0 ),1,0 )
statelevel.disagg$state.cong.prop.sq <- tapply(congruence.data.policy$congruent, congruence.data.policy$state, mean)

#separate states with and without policy change --then present sorted probability of congruence
pdf("Figure_B3.pdf", height=8, width=8)

x.axis.text <- 1.1
y.axis.text <- .7
label.text <- 1.4
title.text <- 1.6
state.text <- .8

prop.cong.no.change.sorted <- sort(statelevel.disagg$state.cong.prop.sq[statelevel.disagg$any.policy.change==0 ])

layout(c(1,2), heights=c(.8,1))
par(mar=c(4,1,1,1))

prop.cong.yes.change.sorted <- sort(statelevel.disagg$state.cong.prop.sq[statelevel.disagg$any.policy.change==1 ])
y <- c(length(prop.cong.yes.change.sorted):1)
plot(prop.cong.yes.change.sorted, y, axes=F, xlab="",ylab="", main="", pch=19, type="n")
text(prop.cong.yes.change.sorted, y, names(prop.cong.yes.change.sorted), cex=state.text)
axis(1, at=seq(0,1,.25), mgp=c(2,.5,0), cex.axis=x.axis.text)
#axis(2, at=y, labels=names(prop.cong.yes.change.sorted), las=1,mgp=c(2,.5,0), cex.axis=y.axis.text)
mtext("Probability that policy was congruent with opinion",1, line=2, cex=label.text)
mtext("States with policy change", 3, line=-2,adj=.9, cex=title.text,font=2)
box()

y <- c(length(prop.cong.no.change.sorted):1)
plot(prop.cong.no.change.sorted, y, axes=F, xlab="",ylab="", main="", pch=19, type="n")
text(prop.cong.no.change.sorted, y, names(prop.cong.no.change.sorted), cex=state.text)
axis(1, at=seq(0,1,.25), mgp=c(2,.5,0), cex.axis=x.axis.text)
#axis(2, at=y, labels=names(prop.cong.no.change.sorted), las=1,mgp=c(2,.5,0), cex.axis=y.axis.text)
mtext("Probability that policy was congruent with opinion",1, line=2, cex=label.text)
mtext("States with no policy change", 3, line=-3,cex=title.text,font=2)
box()

dev.off()
#CREATE 2x2 Table at Bottom of Figure B-3
#another way -- do 2*2 table with uncertainty (as in KLMP). In each sim, get each cell. then build vec
no.policy.change.disagg.opinion.no.policy.change.disagg <- rep(NA, n.sims)
yes.policy.change.disagg.opinion.no.policy.change.disagg <- rep(NA, n.sims)
no.policy.change.disagg.opinion.yes.policy.change.disagg <- rep(NA, n.sims)
yes.policy.change.disagg.opinion.yes.policy.change.disagg <- rep(NA, n.sims)

for (i in 1:n.sims){
	ok <- congruence.data.policy$sim ==paste0("X", i)
 	foo <-	table( congruence.data.policy$any.policy.change[ok] , congruence.data.policy$sq.sims[ok] >= .5)
	no.policy.change.disagg.opinion.no.policy.change.disagg[i] <- foo[1] / length(congruence.data.policy$any.policy.change[ok])
	yes.policy.change.disagg.opinion.no.policy.change.disagg[i] <- foo[2] / length(congruence.data.policy$any.policy.change[ok])
	no.policy.change.disagg.opinion.yes.policy.change.disagg[i] <- foo[3] / length(congruence.data.policy$any.policy.change[ok])
	yes.policy.change.disagg.opinion.yes.policy.change.disagg[i] <- foo[4] / length(congruence.data.policy$any.policy.change[ok])
}
quantile(	no.policy.change.disagg.opinion.no.policy.change.disagg, c(.025,.5,.975))
quantile(	yes.policy.change.disagg.opinion.no.policy.change.disagg, c(.025,.5,.975))
quantile(	no.policy.change.disagg.opinion.yes.policy.change.disagg, c(.025,.5,.975))
quantile(	yes.policy.change.disagg.opinion.yes.policy.change.disagg, c(.025,.5,.975))


########################################################################################
#now bring in courts
########################################################################################
#MRP
##########################################################
#first, focusing  on reforming the status quo at all.

#plot the probability of a court decision overturning a law as a function of public opinion
#code states that didnt' hear cases as zero
table(statelevel.mrp$abortion.policy.change, statelevel.mrp$strike.overall)# CA AND FL
table(statelevel.mrp$state[statelevel.mrp$abortion.policy.change==1 & 	statelevel.mrp$strike.overall==1])
table(statelevel.mrp$state[statelevel.mrp$abortion.policy.change==1 & statelevel.mrp$strike.any==1])
#first, which states heard opinions (note in some case judicial action was cutoff by legislative action, as in NY)
model.heard.case.mrp <- glm(heard.case~
		sq.mrp.pred, family=binomial(link="logit"), data=statelevel.mrp)
summary(model.heard.case.mrp)#no relationship

#choose which opinion
x <- statelevel.mrp$sq.mrp.pred

#plot pr(repeal) against public opinion, in all states 
#NOTE IN TEXT NO RELATIONSHIP BETWEEN PUBLIC OPINION AND STATES IN WHICH CASES WERE HEARD
#1st, entire law
model.invalidate.all.states.mrp <- glm(strike.overall ~
		x, family=binomial(link="logit"), data=statelevel.mrp)
summary(model.invalidate.all.states.mrp,detail=T)
#try just among those states that did not pass any policy change
model.invalidate.no.policy.change.mrp <- glm(strike.overall ~
		x, family=binomial(link="logit"), data=statelevel.mrp,
		subset=abortion.policy.change == 0)
summary(model.invalidate.no.policy.change.mrp, detail=T)

#just in states where cases were heard.
model.invalidate.just.heard.mrp <- glm(strike.overall ~
		x, family=binomial(link="logit"), data=statelevel.mrp, subset=heard.case == 1)
summary(model.invalidate.just.heard.mrp,detail=T)

#just look in federal courts
model.invalidate.just.fed.mrp <- glm(strike.overall ~
		x, family=binomial(link="logit"), data=statelevel.mrp, subset=fed.case == 1)
summary(model.invalidate.just.fed.mrp,detail=T)

#ADD UNCERTAINTY HERE
#DO ALL SIMS HERE TO SAVE TIME
courts.all.states.mrp.coefs <- matrix(NA, nrow=n.sims, ncol=2)#do intercept, then slope
courts.no.policy.change.mrp.coefs <- matrix(NA, nrow=n.sims, ncol=2)#do intercept, then slope
courts.states.heard.mrp.coefs <- matrix(NA, nrow=n.sims, ncol=2)#do intercept, then slope
#all.states.intercepts <- rep(NA, n.sims)
for (i in 1:n.sims){
 temp.x <- which(names(statelevel.mrp)=="X1")
 foo <-	glm(statelevel.mrp$strike.overall~statelevel.mrp[,temp.x+i-1],family=binomial(link="logit"))#L#LOOP OVER COLUMNNS WITH SIMMED OPINION
 foo2 <- sim(foo, n.sims=1)
 courts.all.states.mrp.coefs[i,1] <- coef(foo2)[1]#intercept
 courts.all.states.mrp.coefs[i,2] <- coef(foo2)[2]#slope
 ok <- statelevel.mrp$any.policy.change==0
 foo <-	glm(statelevel.mrp$strike.overall[ok]~statelevel.mrp[,temp.x+i-1][ok],family=binomial(link="logit"))#L#LOOP OVER COLUMNNS WITH SIMMED OPINION
 foo2 <- sim(foo, n.sims=1)
 courts.no.policy.change.mrp.coefs[i,1] <- coef(foo2)[1]#intercept
 courts.no.policy.change.mrp.coefs[i,2] <- coef(foo2)[2]#slope
 ok <- statelevel.mrp$heard.case==1
 foo <-	glm(statelevel.mrp$strike.overall[ok]~statelevel.mrp[,temp.x+i-1][ok],family=binomial(link="logit"))#L#LOOP OVER COLUMNNS WITH SIMMED OPINION
 foo2 <- sim(foo, n.sims=1)
 courts.states.heard.mrp.coefs[i,1] <- coef(foo2)[1]#intercept
 courts.states.heard.mrp.coefs[i,2] <- coef(foo2)[2]#slope
}
quantile(courts.all.states.mrp.coefs[,2],c(.025,.5,.975))
quantile(courts.no.policy.change.mrp.coefs[,2],c(.05,.5,.975))
quantile(courts.states.heard.mrp.coefs[,2],c(.05,.5,.975))
#what is the signficiance of each slope curve?
mean(courts.all.states.mrp.coefs[,2]>0)
mean(courts.no.policy.change.mrp.coefs[,2]>0)
mean(courts.states.heard.mrp.coefs[,2]>0)

#get 95% confidence interval estimates
all.states.lower.mrp.coefs <- c( quantile(courts.all.states.mrp.coefs[,1],(.025)) ,   quantile(courts.all.states.mrp.coefs[,2],(.025)) )

#NOW MAKE Plot 
pdf("Figure_6.pdf", height=9,width=7)

#par(mfrow=c(3,1))
layout(c(1,2,3,4), heights=c(1,1,1,.1))
par(mar=c(2,5,.5,2))
axis.text <- 1.2
title.text <- 1.5
label.text <- 1.3
model.text <- 1.2
states.text <- 1
state.label.adjust <- -.05
form.x <- .7#coordinates for logit equation
form.y <- .3
n.lines <- 1000 #number of simulations to show in plots
offset <- .03

#1)invalidate entire of all states
y <-  statelevel.mrp$strike.overall
x <- statelevel.mrp$sq.mrp.pred
ok <- T
plot(x[ok],y[ok], type="n", ylim=c(0,1.05), xlim=c(.2,.8),axes=F, xlab="", ylab="")
axis(1,  at=seq(0,1,.10), labels=seq(0,100,10), mgp=c(2,.5,0), cex.axis=axis.text, las=1)
axis(2, mgp=c(2,.5,0), cex.axis=axis.text, las=1)

#get 95% confidence intervals at every level of opinion
opinion.sequence <- seq(.2,.8, .05)
courts.all.states.preds.mrp <- matrix(NA, nrow=length(opinion.sequence), ncol=3)
for (j in 1:length(opinion.sequence)){
	courts.all.states.preds.mrp[j,1] <-  quantile(invlogit(courts.all.states.mrp.coefs[,1] + courts.all.states.mrp.coefs[,2]*opinion.sequence[j]), .025)
	courts.all.states.preds.mrp[j,2]  <- quantile(invlogit(courts.all.states.mrp.coefs[,1] + courts.all.states.mrp.coefs[,2]*opinion.sequence[j]), .5)
	courts.all.states.preds.mrp[j,3]  <- quantile(invlogit(courts.all.states.mrp.coefs[,1] + courts.all.states.mrp.coefs[,2]*opinion.sequence[j]), .975)
}
polygon(c(rev(opinion.sequence), opinion.sequence), c(rev(courts.all.states.preds.mrp[ ,3]), courts.all.states.preds.mrp[ ,1]), col = "gray70", border = NA)
points(opinion.sequence, courts.all.states.preds.mrp[ ,2], col="black", type="l")

#plot every other state (in order of x) up and down
#states with invalidations
ok2 <- y==1
x.sorted <- sort(x[ok2]) 
state.sorted <- statelevel.mrp$state[ok2][order(x[ok2])]
temp <- seq(1, length(state.sorted), 1)
odd <- which (temp %% 2 == 1 )
state.sorted.1 <- state.sorted[odd]
even <- which (temp %% 2 != 1 )
state.sorted.2 <- state.sorted[even]

text(x.sorted[odd], 1+offset,  state.sorted.1, cex=1, xpd=NA,srt=0)
text(x.sorted[even],1-offset,  state.sorted.2, cex=1, xpd=NA,srt=0)
text(x[y==0], y[y==0],"|", cex=.8)
mtext("A) All states", 3, cex=title.text, adj=.1, line=-2,font=2)#mtext("State-level support for changing status quo", 1, cex=label.text, line=2)
mtext("Pr(judicial invalidation)",2, srt=90, line=3,  cex=label.text)
box()

#2 ) invalidate entire, of states with no policy change
plot(x,y, type="n",ylim=c(0,1.05), xlim=c(.2,.8),axes=F, xlab="", ylab="")
axis(1,  at=seq(0,1,.10), labels=seq(0,100,10), mgp=c(2,.5,0), cex.axis=axis.text, las=1)
axis(2, mgp=c(2,.5,0), cex.axis=axis.text, las=1)

courts.no.policy.change.preds.mrp <- matrix(NA, nrow=length(opinion.sequence), ncol=3)
for (j in 1:length(opinion.sequence)){
	courts.no.policy.change.preds.mrp[j,1] <-  quantile(invlogit(courts.no.policy.change.mrp.coefs[,1] + courts.no.policy.change.mrp.coefs[,2]*opinion.sequence[j]), .025)
	courts.no.policy.change.preds.mrp[j,2]  <- quantile(invlogit(courts.no.policy.change.mrp.coefs[,1] + courts.no.policy.change.mrp.coefs[,2]*opinion.sequence[j]), .5)
	courts.no.policy.change.preds.mrp[j,3]  <- quantile(invlogit(courts.no.policy.change.mrp.coefs[,1] + courts.no.policy.change.mrp.coefs[,2]*opinion.sequence[j]), .975)
}

polygon(c(rev(opinion.sequence), opinion.sequence), c(rev(courts.no.policy.change.preds.mrp[ ,3]), courts.no.policy.change.preds.mrp[ ,1]), col = "gray70", border = NA)
points(opinion.sequence, courts.no.policy.change.preds.mrp[ ,2], col="black", type="l")

#add text
ok <- statelevel.mrp$any.policy.change==0 & y==0
text(x[ok], y[ok],"|", cex=.8)
x.sorted <- sort(x[ok2]) 
state.sorted <- statelevel.mrp$state[ok2][order(x[ok2])]
temp <- seq(1, length(state.sorted), 1)
odd <- which (temp %% 2 == 1 )
state.sorted.1 <- state.sorted[odd]
even <- which (temp %% 2 != 1 )
state.sorted.2 <- state.sorted[even]
text(x.sorted[odd], 1+offset,  state.sorted.1, cex=1, xpd=NA,srt=0)
text(x.sorted[even],1-offset,  state.sorted.2, cex=1, xpd=NA,srt=0)

mtext("B) States with\nno legislative\npolicy change", 3, cex=title.text,  adj=.1, line=-7,font=2)
mtext("Pr(judicial invalidation)",2, srt=90, line=3, cex=label.text)
box()

#3) Among states in which cases were brought
plot(x,y, type="n",ylim=c(0,1.05), xlim=c(.2,.8),axes=F, xlab="", ylab="")
axis(1,  at=seq(0,1,.10), labels=seq(0,100,10), mgp=c(2,.5,0), cex.axis=axis.text, las=1)
axis(2, mgp=c(2,.5,0), cex.axis=axis.text, las=1)
courts.states.heard.preds.mrp <- matrix(NA, nrow=length(opinion.sequence), ncol=3)
for (j in 1:length(opinion.sequence)){
	courts.states.heard.preds.mrp[j,1] <-  quantile(invlogit(courts.states.heard.mrp.coefs[,1] + courts.states.heard.mrp.coefs[,2]*opinion.sequence[j]), .025)
	courts.states.heard.preds.mrp[j,2]  <- quantile(invlogit(courts.states.heard.mrp.coefs[,1] + courts.states.heard.mrp.coefs[,2]*opinion.sequence[j]), .5)
	courts.states.heard.preds.mrp[j,3]  <- quantile(invlogit(courts.states.heard.mrp.coefs[,1] + courts.states.heard.mrp.coefs[,2]*opinion.sequence[j]), .975)
}

polygon(c(rev(opinion.sequence), opinion.sequence), c(rev(courts.states.heard.preds.mrp[ ,3]), courts.states.heard.preds.mrp[ ,1]), col = "gray70", border = NA)
points(opinion.sequence, courts.states.heard.preds.mrp[ ,2], col="black", type="l")

ok <- statelevel.mrp$heard.case==1 & y==0
text(x[ok], y[ok],"|", cex=1)
ok2 <-  statelevel.mrp$heard.case==1 & y==1
x.sorted <- sort(x[ok2]) 
state.sorted <- statelevel.mrp$state[ok2][order(x[ok2])]
temp <- seq(1, length(state.sorted), 1)
odd <- which (temp %% 2 == 1 )
state.sorted.1 <- state.sorted[odd]
even <- which (temp %% 2 != 1 )
state.sorted.2 <- state.sorted[even]
text(x.sorted[odd], 1+offset,  state.sorted.1, cex=1, xpd=NA,srt=0)
text(x.sorted[even],1-offset,  state.sorted.2, cex=1, xpd=NA,srt=0)

mtext("C) States where\ncases were heard", 3, cex=title.text, adj=.1, line=-5,font=2)
mtext("Pr(judicial invalidation)",2, srt=90, line=3, cex=label.text)
box()

par(mar=c(0,0,0,0))
plot(0,0, type="n", xlab="", ylab="", main="", axes=F)
mtext("State-level support for changing status quo", 1, cex=label.text, line=-2)

dev.off()

#what is the substantive effect of moving from 40\% to 60\% opinion, based on model looking only at states where no policy change occcurre.d
invlogit(median(courts.no.policy.change.mrp.coefs[,1]) + median(courts.no.policy.change.mrp.coefs[,2])*.4)
invlogit(median(courts.no.policy.change.mrp.coefs[,1]) + median(courts.no.policy.change.mrp.coefs[,2])*.6)

#NOW WITH DISAGG
#first, focusing  on reforming the status quo at all.

#plot the probability of a court decision overturning a law as a function of public opinion
#code states that didnt' hear cases as zero
table(statelevel.disagg$abortion.policy.change, statelevel.disagg$strike.overall)# CA AND FL
table(statelevel.disagg$state[statelevel.disagg$abortion.policy.change==1 & 	statelevel.disagg$strike.overall==1])
table(statelevel.disagg$state[statelevel.disagg$abortion.policy.change==1 & statelevel.disagg$strike.any==1])
#first, which states heard opinions (note in some case judicial action was cutoff by legislative action, as in NY)
model.heard.case.disagg <- glm(heard.case~
		sq.disagg.pred, family=binomial(link="logit"), data=statelevel.disagg)
summary(model.heard.case.disagg)#no relationship

#choose which opinion
x <- statelevel.disagg$sq.disagg.pred

#plot pr(repeal) against public opinion, in all states 
#NOTE IN TEXT NO RELATIONSHIP BETWEEN PUBLIC OPINION AND STATES IN WHICH CASES WERE HEARD
#1st, entire law
model.invalidate.all.states.disagg <- glm(strike.overall ~
		x, family=binomial(link="logit"), data=statelevel.disagg)
summary(model.invalidate.all.states.disagg,detail=T)
#try just among those states that did not pass any policy change
model.invalidate.no.policy.change.disagg <- glm(strike.overall ~
		x, family=binomial(link="logit"), data=statelevel.disagg,
		subset=abortion.policy.change == 0)
summary(model.invalidate.no.policy.change.disagg, detail=T)

#just in states where cases were heard.
model.invalidate.just.heard.disagg <- glm(strike.overall ~
		x, family=binomial(link="logit"), data=statelevel.disagg, subset=heard.case == 1)
summary(model.invalidate.just.heard.disagg,detail=T)

#just look in federal courts
model.invalidate.just.fed.disagg <- glm(strike.overall ~
		x, family=binomial(link="logit"), data=statelevel.disagg, subset=fed.case == 1)
summary(model.invalidate.just.fed.disagg,detail=T)

#ADD UNCERTAINTY HERE
#DO ALL SIMS HERE TO SAVE TIME
courts.all.states.disagg.coefs <- matrix(NA, nrow=n.sims, ncol=2)#do intercept, then slope
courts.no.policy.change.disagg.coefs <- matrix(NA, nrow=n.sims, ncol=2)#do intercept, then slope
courts.states.heard.disagg.coefs <- matrix(NA, nrow=n.sims, ncol=2)#do intercept, then slope
#all.states.intercepts <- rep(NA, n.sims)
for (i in 1:n.sims){
 temp.x <- which(names(statelevel.disagg)=="X1")
 foo <-	glm(statelevel.disagg$strike.overall~statelevel.disagg[,temp.x+i-1],family=binomial(link="logit"))#L#LOOP OVER COLUMNNS WITH SIMMED OPINION
 foo2 <- sim(foo, n.sims=1)
 courts.all.states.disagg.coefs[i,1] <- coef(foo2)[1]#intercept
 courts.all.states.disagg.coefs[i,2] <- coef(foo2)[2]#slope
 ok <- statelevel.disagg$any.policy.change==0
 foo <-	glm(statelevel.disagg$strike.overall[ok]~statelevel.disagg[,temp.x+i-1][ok],family=binomial(link="logit"))#L#LOOP OVER COLUMNNS WITH SIMMED OPINION
 foo2 <- sim(foo, n.sims=1)
 courts.no.policy.change.disagg.coefs[i,1] <- coef(foo2)[1]#intercept
 courts.no.policy.change.disagg.coefs[i,2] <- coef(foo2)[2]#slope
 ok <- statelevel.disagg$heard.case==1
 foo <-	glm(statelevel.disagg$strike.overall[ok]~statelevel.disagg[,temp.x+i-1][ok],family=binomial(link="logit"))#L#LOOP OVER COLUMNNS WITH SIMMED OPINION
 foo2 <- sim(foo, n.sims=1)
 courts.states.heard.disagg.coefs[i,1] <- coef(foo2)[1]#intercept
 courts.states.heard.disagg.coefs[i,2] <- coef(foo2)[2]#slope
}
quantile(courts.all.states.disagg.coefs[,2],c(.025,.5,.975))
quantile(courts.no.policy.change.disagg.coefs[,2],c(.05,.5,.975))
quantile(courts.states.heard.disagg.coefs[,2],c(.05,.5,.975))
mean(courts.all.states.disagg.coefs[,2]>0)
mean(courts.no.policy.change.disagg.coefs[,2]>0)
mean(courts.states.heard.disagg.coefs[,2]>0)

#get 95% confidence interval estimates
all.states.lower.disagg.coefs <- c( quantile(courts.all.states.disagg.coefs[,1],(.025)) ,   quantile(courts.all.states.disagg.coefs[,2],(.025)) )

#NOW MAKE Plot 
pdf("Figure_B4.pdf", height=9,width=7)

#par(mfrow=c(3,1))
layout(c(1,2,3,4), heights=c(1,1,1,.1))
par(mar=c(2,5,.5,2))
axis.text <- 1.2
title.text <- 1.5
label.text <- 1.3
model.text <- 1.2
states.text <- 1
state.label.adjust <- -.05
form.x <- .7#coordinates for logit equation
form.y <- .3
n.lines <- 1000 #number of simulations to show in plots
offset <- .03

#1)invalidate entire of all states
y <-  statelevel.disagg$strike.overall
x <- statelevel.disagg$sq.disagg.pred
ok <- T
plot(x[ok],y[ok], type="n", ylim=c(0,1.05), xlim=c(.2,.8),axes=F, xlab="", ylab="")
axis(1,  at=seq(0,1,.10), labels=seq(0,100,10), mgp=c(2,.5,0), cex.axis=axis.text, las=1)
axis(2, mgp=c(2,.5,0), cex.axis=axis.text, las=1)

#get 95% confidence intervals at every level of opinion
opinion.sequence <- seq(.2,.8, .05)
courts.all.states.preds.disagg <- matrix(NA, nrow=length(opinion.sequence), ncol=3)
for (j in 1:length(opinion.sequence)){
	courts.all.states.preds.disagg[j,1] <-  quantile(invlogit(courts.all.states.disagg.coefs[,1] + courts.all.states.disagg.coefs[,2]*opinion.sequence[j]), .025)
	courts.all.states.preds.disagg[j,2]  <- quantile(invlogit(courts.all.states.disagg.coefs[,1] + courts.all.states.disagg.coefs[,2]*opinion.sequence[j]), .5)
	courts.all.states.preds.disagg[j,3]  <- quantile(invlogit(courts.all.states.disagg.coefs[,1] + courts.all.states.disagg.coefs[,2]*opinion.sequence[j]), .975)
}
polygon(c(rev(opinion.sequence), opinion.sequence), c(rev(courts.all.states.preds.disagg[ ,3]), courts.all.states.preds.disagg[ ,1]), col = "gray70", border = NA)
points(opinion.sequence, courts.all.states.preds.disagg[ ,2], col="black", type="l")

#plot every other state (in order of x) up and down
#states with invalidations
ok2 <- y==1
x.sorted <- sort(x[ok2]) 
state.sorted <- statelevel.disagg$state[ok2][order(x[ok2])]
temp <- seq(1, length(state.sorted), 1)
odd <- which (temp %% 2 == 1 )
state.sorted.1 <- state.sorted[odd]
even <- which (temp %% 2 != 1 )
state.sorted.2 <- state.sorted[even]

text(x.sorted[odd], 1+offset,  state.sorted.1, cex=1, xpd=NA,srt=0)
text(x.sorted[even],1-offset,  state.sorted.2, cex=1, xpd=NA,srt=0)
text(x[y==0], y[y==0],"|", cex=.8)
mtext("A) All states", 3, cex=title.text, adj=.1, line=-2,font=2)
#mtext("State-level support for changing status quo", 1, cex=label.text, line=2)
mtext("Pr(judicial invalidation)",2, srt=90, line=3,  cex=label.text)
box()

#2 ) invalidate entire, of states with no policy change
plot(x,y, type="n",ylim=c(0,1.05), xlim=c(.2,.8),axes=F, xlab="", ylab="")
axis(1,  at=seq(0,1,.10), labels=seq(0,100,10), mgp=c(2,.5,0), cex.axis=axis.text, las=1)
axis(2, mgp=c(2,.5,0), cex.axis=axis.text, las=1)

courts.no.policy.change.preds.disagg <- matrix(NA, nrow=length(opinion.sequence), ncol=3)
for (j in 1:length(opinion.sequence)){
	courts.no.policy.change.preds.disagg[j,1] <-  quantile(invlogit(courts.no.policy.change.disagg.coefs[,1] + courts.no.policy.change.disagg.coefs[,2]*opinion.sequence[j]), .025)
	courts.no.policy.change.preds.disagg[j,2]  <- quantile(invlogit(courts.no.policy.change.disagg.coefs[,1] + courts.no.policy.change.disagg.coefs[,2]*opinion.sequence[j]), .5)
	courts.no.policy.change.preds.disagg[j,3]  <- quantile(invlogit(courts.no.policy.change.disagg.coefs[,1] + courts.no.policy.change.disagg.coefs[,2]*opinion.sequence[j]), .975)
}

polygon(c(rev(opinion.sequence), opinion.sequence), c(rev(courts.no.policy.change.preds.disagg[ ,3]), courts.no.policy.change.preds.disagg[ ,1]), col = "gray70", border = NA)
points(opinion.sequence, courts.no.policy.change.preds.disagg[ ,2], col="black", type="l")

#add text
ok <- statelevel.disagg$any.policy.change==0 & y==0
text(x[ok], y[ok],"|", cex=.8)
x.sorted <- sort(x[ok2]) 
state.sorted <- statelevel.disagg$state[ok2][order(x[ok2])]
temp <- seq(1, length(state.sorted), 1)
odd <- which (temp %% 2 == 1 )
state.sorted.1 <- state.sorted[odd]
even <- which (temp %% 2 != 1 )
state.sorted.2 <- state.sorted[even]
text(x.sorted[odd], 1+offset,  state.sorted.1, cex=1, xpd=NA,srt=0)
text(x.sorted[even],1-offset,  state.sorted.2, cex=1, xpd=NA,srt=0)

mtext("B) States with\nno legislative\npolicy change", 3, cex=title.text,  adj=.1, line=-7,font=2)
#mtext("State-level support for changing status quo", 1, cex=label.text, line=2)
mtext("Pr(judicial invalidation)",2, srt=90, line=3, cex=label.text)
box()

#3) Among states in which cases were brought
plot(x,y, type="n",ylim=c(0,1.05), xlim=c(.2,.8),axes=F, xlab="", ylab="")
axis(1,  at=seq(0,1,.10), labels=seq(0,100,10), mgp=c(2,.5,0), cex.axis=axis.text, las=1)
axis(2, mgp=c(2,.5,0), cex.axis=axis.text, las=1)
courts.states.heard.preds.disagg <- matrix(NA, nrow=length(opinion.sequence), ncol=3)
for (j in 1:length(opinion.sequence)){
	courts.states.heard.preds.disagg[j,1] <-  quantile(invlogit(courts.states.heard.disagg.coefs[,1] + courts.states.heard.disagg.coefs[,2]*opinion.sequence[j]), .025)
	courts.states.heard.preds.disagg[j,2]  <- quantile(invlogit(courts.states.heard.disagg.coefs[,1] + courts.states.heard.disagg.coefs[,2]*opinion.sequence[j]), .5)
	courts.states.heard.preds.disagg[j,3]  <- quantile(invlogit(courts.states.heard.disagg.coefs[,1] + courts.states.heard.disagg.coefs[,2]*opinion.sequence[j]), .975)
}

polygon(c(rev(opinion.sequence), opinion.sequence), c(rev(courts.states.heard.preds.disagg[ ,3]), courts.states.heard.preds.disagg[ ,1]), col = "gray70", border = NA)
points(opinion.sequence, courts.states.heard.preds.disagg[ ,2], col="black", type="l")

ok <- statelevel.disagg$heard.case==1 & y==0
text(x[ok], y[ok],"|", cex=1)
ok2 <-  statelevel.disagg$heard.case==1 & y==1
x.sorted <- sort(x[ok2]) 
state.sorted <- statelevel.disagg$state[ok2][order(x[ok2])]
temp <- seq(1, length(state.sorted), 1)
odd <- which (temp %% 2 == 1 )
state.sorted.1 <- state.sorted[odd]
even <- which (temp %% 2 != 1 )
state.sorted.2 <- state.sorted[even]
text(x.sorted[odd], 1+offset,  state.sorted.1, cex=1, xpd=NA,srt=0)
text(x.sorted[even],1-offset,  state.sorted.2, cex=1, xpd=NA,srt=0)

mtext("C) States where\ncases were heard", 3, cex=title.text, adj=.1, line=-5,font=2)
#mtext("State-level support for changing status quo", 1, cex=label.text, line=2)
mtext("Pr(judicial invalidation)",2, srt=90, line=3, cex=label.text)
box()

par(mar=c(0,0,0,0))
plot(0,0, type="n", xlab="", ylab="", main="", axes=F)
mtext("State-level support for changing status quo", 1, cex=label.text, line=-2)

dev.off()



#NOW DO CONGRUENCE ANALYSIS, focusing only on states where cases were heard
#MRP
ok <- statelevel.mrp$heard.case==1# & statelevel.mrp$strike.overall==1
table(statelevel.mrp$strike.overall[ok], statelevel.mrp$sq.mrp.pred[ok]>.5)
table(statelevel.mrp$strike.overall, statelevel.mrp$sq.mrp.pred>.5)

data.frame(statelevel.mrp$state[ok], statelevel.mrp$strike.overall[ok], statelevel.mrp$sq.mrp.pred[ok])[order(statelevel.mrp$strike.overall[ok],statelevel.mrp$sq.mrp.pred[ok]),]

#just keep state, any policy change and sims
temp.end <- paste0("X",n.sims)
temp1 <- statelevel.mrp[,which(names(statelevel.mrp)=="X1"): which(names(statelevel.mrp)==paste0("X",n.sims)) ]
strike.overall <- statelevel.mrp[,"strike.overall"]
state <- statelevel.mrp[,"state"]
heard.case <- statelevel.mrp[,"heard.case"]
temp.df <- data.frame(cbind(state, strike.overall ,heard.case, temp1))
temp.df <- subset(temp.df,heard.case==1)
temp.df$heard.case<-NULL
temp.df$state <- droplevels(temp.df$state)#drop missing levels

congruence.data.courts <- melt(temp.df, id.vars = c("state", "strike.overall"), value.name="sq.sims", variable.name="sim")
congruence.data.courts <- congruence.data.courts[order(congruence.data.courts$state),]

#now create congruence variable
congruence.data.courts$congruent <- ifelse( ( congruence.data.courts$sq.sims >=.5 & congruence.data.courts$strike.overall == 1 ) |
	 (congruence.data.courts$sq.sims <.5 & congruence.data.courts$strike.overall == 0 ),1,0 )
temp.df$courts.cong.prop.sq <- tapply(congruence.data.courts$congruent, congruence.data.courts$state, mean)

#separate states with and without invalidation --then present sorted probability of congruence
pdf("Figure_7.pdf", height=8, width=8)

x.axis.text <- 1.1
y.axis.text <- .7
label.text <- 1.4
title.text <- 1.6
state.text <- .7

courts.cong.prop.sq.sorted.strike.no <- sort(temp.df$courts.cong.prop.sq [temp.df$strike.overall==0])
y <- c(length(courts.cong.prop.sq.sorted.strike.no):1)

layout(c(1,2), heights=c(1,.8))
par(mar=c(4,1,0,0))
plot(courts.cong.prop.sq.sorted.strike.no, y, axes=F, xlab="",ylab="", main="", pch=19, xlim=c(0,1), type="n")
text(courts.cong.prop.sq.sorted.strike.no, y, names(courts.cong.prop.sq.sorted.strike.no), cex=state.text)
axis(1, at=seq(0,1,.25), mgp=c(2,.5,0), cex.axis=x.axis.text)
#axis(2, at=y, labels=names(courts.cong.prop.sq.sorted.strike.no), las=1,mgp=c(2,.5,0), cex.axis=y.axis.text)
mtext("Probability that court decision was congruent with opinion",1, line=2, cex=label.text)
mtext("States where courts upheld statute", 3, line=-2,cex=title.text,font=2)
courts.cong.prop.sq.sorted.strike.yes <- sort(temp.df$courts.cong.prop.sq [temp.df$strike.overall==1])
y <- c(length(courts.cong.prop.sq.sorted.strike.yes):1)
box()

plot(courts.cong.prop.sq.sorted.strike.yes, y, axes=F, xlab="",ylab="", main="", pch=19,xlim=c(0,1), type="n")
axis(1, at=seq(0,1,.25), mgp=c(2,.5,0), cex.axis=x.axis.text)
text(courts.cong.prop.sq.sorted.strike.yes, y, names(courts.cong.prop.sq.sorted.strike.yes), cex=state.text)
#axis(2, at=y, labels=names(courts.cong.prop.sq.sorted.strike.yes), las=1,mgp=c(2,.5,0), cex.axis=y.axis.text)
mtext("Probability that court decision was congruent with opinion",1, line=2, cex=label.text)
mtext("States where courts invalidated statute", 1, adj=.5, line=-2,cex=title.text,font=2)
box()

dev.off()

#CREATE 2x2 table on bottom of Figure 7
#another way -- do 2*2 table with uncertainty (as in KLMP). In each sim, get each cell. then build vec
courts.uphold.opinion.no.policy.change.mrp <- rep(NA, n.sims)
courts.strike.opinion.no.policy.change.mrp <- rep(NA, n.sims)
courts.uphold.opinion.yes.policy.change.mrp <- rep(NA, n.sims)
courts.strike.opinion.yes.policy.change.mrp <- rep(NA, n.sims)

for (i in 1:n.sims){
	ok <- congruence.data.courts$sim ==paste0("X", i)
 	foo <-	table( congruence.data.courts$strike.overall[ok] , congruence.data.courts$sq.sims[ok] >= .5)
	courts.uphold.opinion.no.policy.change.mrp[i] <- foo[1] / length(congruence.data.courts$strike.overall[ok])
	courts.strike.opinion.no.policy.change.mrp[i] <- foo[2] / length(congruence.data.courts$strike.overall[ok])
	courts.uphold.opinion.yes.policy.change.mrp[i] <- foo[3] / length(congruence.data.courts$strike.overall[ok])
	courts.strike.opinion.yes.policy.change.mrp[i] <- foo[4] / length(congruence.data.courts$strike.overall[ok])
}
quantile(	courts.uphold.opinion.no.policy.change.mrp, c(.025,.5,.975))
quantile(	courts.strike.opinion.no.policy.change.mrp, c(.025,.5,.975))
quantile(	courts.uphold.opinion.yes.policy.change.mrp, c(.025,.5,.975))
quantile(	courts.strike.opinion.yes.policy.change.mrp, c(.025,.5,.975))

####################################################
#DISAGG
########################################################################################################
ok <- statelevel.disagg$heard.case==1# & statelevel.disagg$strike.overall==1
table(statelevel.disagg$strike.overall[ok], statelevel.disagg$sq.disagg.pred[ok]>.5)
table(statelevel.disagg$strike.overall, statelevel.disagg$sq.disagg.pred>.5)

data.frame(statelevel.disagg$state[ok], statelevel.disagg$strike.overall[ok], statelevel.disagg$sq.disagg.pred[ok])[order(statelevel.disagg$strike.overall[ok],statelevel.disagg$sq.disagg.pred[ok]),]

#just keep state, any policy change and sims
temp.end <- paste0("X",n.sims)
temp1 <- statelevel.disagg[,which(names(statelevel.disagg)=="X1"): which(names(statelevel.disagg)==paste0("X",n.sims)) ]
strike.overall <- statelevel.disagg[,"strike.overall"]
state <- statelevel.disagg[,"state"]
heard.case <- statelevel.disagg[,"heard.case"]
temp.df <- data.frame(cbind(state, strike.overall ,heard.case, temp1))
temp.df <- subset(temp.df,heard.case==1)
temp.df$heard.case<-NULL
temp.df$state <- droplevels(temp.df$state)#drop missing levels

congruence.data.courts <- melt(temp.df, id.vars = c("state", "strike.overall"), value.name="sq.sims", variable.name="sim")
congruence.data.courts <- congruence.data.courts[order(congruence.data.courts$state),]

#now create congruence variable
congruence.data.courts$congruent <- ifelse( ( congruence.data.courts$sq.sims >=.5 & congruence.data.courts$strike.overall == 1 ) |
	 (congruence.data.courts$sq.sims <.5 & congruence.data.courts$strike.overall == 0 ),1,0 )
temp.df$courts.cong.prop.sq <- tapply(congruence.data.courts$congruent, congruence.data.courts$state, mean)

#separate states with and without invalidation --then present sorted probability of congruence
pdf("Figure_B5.pdf", height=8, width=8)

x.axis.text <- 1.1
y.axis.text <- .7
label.text <- 1.4
title.text <- 1.6
state.text <- .7

courts.cong.prop.sq.sorted.strike.no <- sort(temp.df$courts.cong.prop.sq [temp.df$strike.overall==0])
y <- c(length(courts.cong.prop.sq.sorted.strike.no):1)

layout(c(1,2), heights=c(1,.8))
par(mar=c(4,1,0,0))
plot(courts.cong.prop.sq.sorted.strike.no, y, axes=F, xlab="",ylab="", main="", pch=19, xlim=c(0,1), type="n")
text(courts.cong.prop.sq.sorted.strike.no, y, names(courts.cong.prop.sq.sorted.strike.no), cex=state.text)
axis(1, at=seq(0,1,.25), mgp=c(2,.5,0), cex.axis=x.axis.text)
#axis(2, at=y, labels=names(courts.cong.prop.sq.sorted.strike.no), las=1,mgp=c(2,.5,0), cex.axis=y.axis.text)
mtext("Probability that court decision was congruent with opinion",1, line=2, cex=label.text)
mtext("States where courts upheld statute", 3, line=-2,cex=title.text,font=2)
courts.cong.prop.sq.sorted.strike.yes <- sort(temp.df$courts.cong.prop.sq [temp.df$strike.overall==1])
y <- c(length(courts.cong.prop.sq.sorted.strike.yes):1)
box()

plot(courts.cong.prop.sq.sorted.strike.yes, y, axes=F, xlab="",ylab="", main="", pch=19,xlim=c(0,1), type="n")
axis(1, at=seq(0,1,.25), mgp=c(2,.5,0), cex.axis=x.axis.text)
text(courts.cong.prop.sq.sorted.strike.yes, y, names(courts.cong.prop.sq.sorted.strike.yes), cex=state.text)
#axis(2, at=y, labels=names(courts.cong.prop.sq.sorted.strike.yes), las=1,mgp=c(2,.5,0), cex.axis=y.axis.text)
mtext("Probability that court decision was congruent with opinion",1, line=2, cex=label.text)
mtext("States where courts invalidated statute", 1, adj=.5, line=-2,cex=title.text,font=2)
box()

dev.off()

#CREATE 2x2 table on bottom of Figure B-5
#another way -- do 2*2 table with uncertainty (as in KLMP). In each sim, get each cell. then build vec
courts.uphold.opinion.no.policy.change.disagg <- rep(NA, n.sims)
courts.strike.opinion.no.policy.change.disagg <- rep(NA, n.sims)
courts.uphold.opinion.yes.policy.change.disagg <- rep(NA, n.sims)
courts.strike.opinion.yes.policy.change.disagg <- rep(NA, n.sims)

for (i in 1:n.sims){
	ok <- congruence.data.courts$sim ==paste0("X", i)
 	foo <-	table( congruence.data.courts$strike.overall[ok] , congruence.data.courts$sq.sims[ok] >= .5)
	courts.uphold.opinion.no.policy.change.disagg[i] <- foo[1] / length(congruence.data.courts$strike.overall[ok])
	courts.strike.opinion.no.policy.change.disagg[i] <- foo[2] / length(congruence.data.courts$strike.overall[ok])
	courts.uphold.opinion.yes.policy.change.disagg[i] <- foo[3] / length(congruence.data.courts$strike.overall[ok])
	courts.strike.opinion.yes.policy.change.disagg[i] <- foo[4] / length(congruence.data.courts$strike.overall[ok])
}
quantile(	courts.uphold.opinion.no.policy.change.disagg, c(.025,.5,.975))
quantile(	courts.strike.opinion.no.policy.change.disagg, c(.025,.5,.975))
quantile(	courts.uphold.opinion.yes.policy.change.disagg, c(.025,.5,.975))
quantile(	courts.strike.opinion.yes.policy.change.disagg, c(.025,.5,.975))








